library(plyr)
library(igraph)
library(fmsb)
library(Hmisc)
library(MASS)
library(stargazer)
library(tidyverse)
library(multcomp)
library(gsynth)
data(gsynth)
library(panelView)

setwd("~/Dropbox/China Huawei Paper/Production Materials/Replication/")
master.data <- read.csv("dataFile_China_Huawei.csv"); head(master.data)

# load broadband and telephone data
library(WDI)
library(countrycode)
wb = WDI(country="all", indicator=c("IT.NET.BBND.P2", "IT.MLT.MAIN.P2"), start=2000,end=2020,extra=T,cache=NULL)
names(wb)[which(names(wb) == "IT.NET.BBND.P2")] <- "broadbandpc" #fixed broadband subscribers (per 100 people)"
names(wb)[which(names(wb) == "IT.MLT.MAIN.P2")] <- "telephonepc" #fixed telephone subscribers (per 100 people)"
wb$cabb = countrycode(as.character(as.vector(wb$iso3c)),origin="iso3c",destination="cowc",warn=T)
wb = wb[wb$cabb %in% master.data$cabb,];nrow(wb) # retain all years, but limit to countries in master.data
wb = wb[c("cabb","year", "broadbandpc","telephonepc")]
nrow(master.data)
master.data = join(master.data,wb,type="full",by=c("cabb","year"))
nrow(master.data) 

# Internet penetration
dat <- master.data
dat$treated = 0
threshold=summary(master.data$internetpen)[3]
dat = as.data.frame(dat %>%
  group_by(cabb) %>%
  dplyr::mutate(treated = ifelse(row_number() == 1, 0, NA),
         treated = ifelse(internetpen > threshold, 1, treated)) %>%
  fill(treated))
dat$treated[is.na(dat$treated)==T] <- 0

sub = dat[dat$v2x_polyarchy>0.42,]
sub = sub[sub$cabb!="MNG",]
sub = sub[sub$cabb!="SSD",]
sub = sub[is.na(sub$v2smgovfilprc)==F & is.na(sub$treated)==F & is.na(sub$v2x_polyarchy)==F & is.na(sub$pres.election)==F & is.na(sub$coup.attempts)==F & is.na(sub$successful.coups)==F,]
sub = sub[is.na(sub$pcGDPcur)==F & is.na(sub$GDPcur)==F & is.na(sub$electricitypc)==F,]
nrow(dat);nrow(sub)

# must drop to min.T0=3
out1 <- gsynth(v2smgovfilprc ~ treated + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc, min.T0=3, data=sub, index=c("cabb","year"), estimator = "mc", se=TRUE, nboots=1000, seed=02139)
plot(out1, type="gap", xlim=c(-5,10), ylim=c(-0.2,1.0), main="Democracies\nSocial Media Filtering")
out1$est.avg
out1$N
out1internetPenetration = out1$est.avg

# must drop to min.T0=3
out2 <- gsynth(v2smgovshut ~ treated + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc, min.T0=3, data=sub, index=c("cabb","year"), estimator = "mc", se=TRUE, nboots=1000, seed=02139)
plot(out2, type="gap", xlim=c(-5,10), ylim=c(-0.2,1.0), main="Democracies\nInternet Shutdowns")
out2$est.avg
out2internetPenetration = out2$est.avg

# must drop to min.T0=3
out3 <- gsynth(v2smgovsmmon ~ treated + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc, min.T0=3, data=sub, index=c("cabb","year"), estimator = "mc", se=TRUE, nboots=1000, seed=02139)
plot(out3, type="gap", xlim=c(-5,10), ylim=c(-0.2,1.0), main="Democracies\nSocial Media Monitoring")
out3$est.avg
out3internetPenetration = out3$est.avg

# must drop to min.T0=3
out4 <- gsynth(v2smarrest ~ treated + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc, min.T0=3, data=sub, index=c("cabb","year"), estimator = "mc", se=TRUE, nboots=1000, seed=02139)
plot(out4, type="gap", xlim=c(-5,10), ylim=c(-0.2,1.0), main="Democracies\nArrests for Online Content")
out4$est.avg
out4internetPenetration = out4$est.avg


# Broadband subscribers 
dat <- master.data
dat$treated = 0
threshold=summary(master.data$broadbandpc)[3]
dat = as.data.frame(dat %>%
  group_by(cabb) %>%
  dplyr::mutate(treated = ifelse(row_number() == 1, 0, NA),
         treated = ifelse(broadbandpc > threshold, 1, treated)) %>%
  fill(treated))
dat$treated[is.na(dat$treated)==T] <- 0

sub = dat[dat$v2x_polyarchy>0.42,]
sub = sub[sub$cabb!="MNG",]
sub = sub[sub$cabb!="SSD",]
sub = sub[is.na(sub$v2smgovfilprc)==F & is.na(sub$treated)==F & is.na(sub$v2x_polyarchy)==F & is.na(sub$pres.election)==F & is.na(sub$coup.attempts)==F & is.na(sub$successful.coups)==F,]
sub = sub[is.na(sub$pcGDPcur)==F & is.na(sub$GDPcur)==F & is.na(sub$electricitypc)==F,]
nrow(dat);nrow(sub)

# must drop to min.T0=3
out1 <- gsynth(v2smgovfilprc ~ treated + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc, min.T0=3, data=sub, index=c("cabb","year"), estimator = "mc", se=TRUE, nboots=1000, seed=02139)
plot(out1, type="gap", xlim=c(-5,10), ylim=c(-0.2,1.0), main="Democracies\nSocial Media Filtering")
out1$est.avg
out1broadbandSubscribers = out1$est.avg

# must drop to min.T0=3
out2 <- gsynth(v2smgovshut ~ treated + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc, min.T0=3, data=sub, index=c("cabb","year"), estimator = "mc", se=TRUE, nboots=1000, seed=02139)
plot(out2, type="gap", xlim=c(-5,10), ylim=c(-0.2,1.0), main="Democracies\nInternet Shutdowns")
out2$est.avg
out2broadbandSubscribers = out2$est.avg

# must drop to min.T0=3
out3 <- gsynth(v2smgovsmmon ~ treated + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc, min.T0=3, data=sub, index=c("cabb","year"), estimator = "mc", se=TRUE, nboots=1000, seed=02139)
plot(out3, type="gap", xlim=c(-5,10), ylim=c(-0.2,1.0), main="Democracies\nSocial Media Monitoring")
out3$est.avg
out3broadbandSubscribers = out3$est.avg

# must drop to min.T0=3
out4 <- gsynth(v2smarrest ~ treated + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc, min.T0=3, data=sub, index=c("cabb","year"), estimator = "mc", se=TRUE, nboots=1000, seed=02139)
plot(out4, type="gap", xlim=c(-5,10), ylim=c(-0.2,1.0), main="Democracies\nArrests for Online Content")
out4$est.avg
out4broadbandSubscribers = out4$est.avg


# Telephone subscribers
dat <- master.data
dat$treated = 0
threshold=summary(master.data$telephonepc)[3]
dat = as.data.frame(dat %>%
  group_by(cabb) %>%
  dplyr::mutate(treated = ifelse(row_number() == 1, 0, NA),
         treated = ifelse(telephonepc > threshold, 1, treated)) %>%
  fill(treated))
dat$treated[is.na(dat$treated)==T] <- 0

sub = dat[dat$v2x_polyarchy>0.42,]
sub = sub[sub$cabb!="MNG",]
sub = sub[sub$cabb!="SSD",]
sub = sub[is.na(sub$v2smgovfilprc)==F & is.na(sub$treated)==F & is.na(sub$v2x_polyarchy)==F & is.na(sub$pres.election)==F & is.na(sub$coup.attempts)==F & is.na(sub$successful.coups)==F,]
sub = sub[is.na(sub$pcGDPcur)==F & is.na(sub$GDPcur)==F & is.na(sub$electricitypc)==F,]
nrow(dat);nrow(sub)

# back to min.T0=5 for these models
out1 <- gsynth(v2smgovfilprc ~ treated + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc, min.T0=5, data=sub, index=c("cabb","year"), estimator = "mc", se=TRUE, nboots=1000, seed=02139)
plot(out1, type="gap", xlim=c(-5,10), ylim=c(-0.2,1.0), main="Democracies\nSocial Media Filtering")
out1$est.avg
out1telephoneSubscribers = out1$est.avg

out2 <- gsynth(v2smgovshut ~ treated + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc, min.T0=5, data=sub, index=c("cabb","year"), estimator = "mc", se=TRUE, nboots=1000, seed=02139)
plot(out2, type="gap", xlim=c(-5,10), ylim=c(-0.2,1.0), main="Democracies\nInternet Shutdowns")
out2$est.avg
out2telephoneSubscribers = out2$est.avg

out3 <- gsynth(v2smgovsmmon ~ treated + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc, min.T0=5, data=sub, index=c("cabb","year"), estimator = "mc", se=TRUE, nboots=1000, seed=02139)
plot(out3, type="gap", xlim=c(-5,10), ylim=c(-0.2,1.0), main="Democracies\nSocial Media Monitoring")
out3$est.avg
out3telephoneSubscribers = out3$est.avg

out4 <- gsynth(v2smarrest ~ treated + v2x_polyarchy + pres.election + coup.attempts + successful.coups + icews.protest + icews.repression + GDPcur + pcGDPcur + electricitypc, min.T0=5, data=sub, index=c("cabb","year"), estimator = "mc", se=TRUE, nboots=1000, seed=02139)
plot(out4, type="gap", xlim=c(-5,10), ylim=c(-0.2,1.0), main="Democracies\nArrests for Online Content")
out4$est.avg
out4telephoneSubscribers = out4$est.avg


# export table
ti = data.frame(rbind(out1internetPenetration,out1broadbandSubscribers,out1telephoneSubscribers))
ti$term = c("Internet Penetration (Value \u2265 Median)","Broadband Subscribers (Value \u2265 Median)","Telephone Subscribers (Value \u2265 Median)")
ti = ti[,c("term","Estimate","S.E.","p.value")]
names(ti) = c("term","estimate","std.error","p.value")
out1 = list(tidy=ti)
class(out1) = "modelsummary_list"

ti = data.frame(rbind(out2internetPenetration,out2broadbandSubscribers,out2telephoneSubscribers))
ti$term = c("Internet Penetration (Value \u2265 Median)","Broadband Subscribers (Value \u2265 Median)","Telephone Subscribers (Value \u2265 Median)")
ti = ti[,c("term","Estimate","S.E.","p.value")]
names(ti) = c("term","estimate","std.error","p.value")
out2 = list(tidy=ti)
class(out2) = "modelsummary_list"

ti = data.frame(rbind(out3internetPenetration,out3broadbandSubscribers,out3telephoneSubscribers))
ti$term = c("Internet Penetration (Value \u2265 Median)","Broadband Subscribers (Value \u2265 Median)","Telephone Subscribers (Value \u2265 Median)")
ti = ti[,c("term","Estimate","S.E.","p.value")]
names(ti) = c("term","estimate","std.error","p.value")
out3 = list(tidy=ti)
class(out3) = "modelsummary_list"

ti = data.frame(rbind(out4internetPenetration,out4broadbandSubscribers,out4telephoneSubscribers))
ti$term = c("Internet Penetration (Value \u2265 Median)","Broadband Subscribers (Value \u2265 Median)","Telephone Subscribers (Value \u2265 Median)")
ti = ti[,c("term","Estimate","S.E.","p.value")]
names(ti) = c("term","estimate","std.error","p.value")
out4 = list(tidy=ti)
class(out4) = "modelsummary_list"


myrows = data.frame(list(c("Country Fixed Effects",rep("Yes",4))
                         ,c("Year Fixed Effects",rep("Yes",4))
                         ,c("Control Variables",rep("Yes",4))
))

export = modelsummary(list("Internet Filtering" = out1
                           , "Internet Shutdowns" = out2
                           ,"Social Media Monitoring"=out3
                           ,"Arrests for Political Content"=out4
)
,add_rows=data.frame(t(myrows))
, stars=c('*' = .1, '**' = 0.05, '***' = .01)
, output="latex_tabular"
)
export = gsub("≥", "$\\\\ge$", export)
export



